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We devise and examine two procrustean entanglement concentration schemes using Quantum Non- 
Demolition (QND) interaction Hamiltonians in the continuous variable regime, applicable for light, 
for atomic ensembles or in a hybrid setting. We thus expand the standard entanglement distillation 
toolbox to the use of a much more general, versatile and experimentally feasible interaction class. 
The first protocol uses Gaussian ancillary modes and a non-Gaussian post-measurement, the second 
a non-Gaussian ancillary mode and a Gaussian post-measurement. We explicitly calculate the 
density matrix elements of the non-Gaussian mixed states resulting from these protocols using an 
elegant Wigner-function based method in a numerically efficient manner. We then quantify the 
entanglement increase calculating the Logarithmic Negativity of the output state and discuss and 
compare the performance of the protocols. 



INTRODUCTION 



One of the biggest hurdles quantum communication 
has yet to overcome is the difficulty of generating highly 
entangled states between distant nodes. Due to the in- 
herent noise in any optical setup, in reality it is only ever 
possible to generate partially entangled states at distant 
locations. The remedy to this problem is entanglement 
distillation [TH3j. Local operations and classical commu- 
nication (LOGO) are used to distil at least one maximally 
entangled bipartite state from multiple weakly entangled 
states. 

One of the latest developments in quantum informa- 
tion science is the use of infinite-dimensional systems or 
continuous variables (CV), such as the amplitude and 
phase quadratures of optical modes or the collective spin 
of macroscopic atomic ensembles. CV systems offer the 
advantage of cheap resources, easy generation and con- 
trol, and off-the-shelf telecommunication components can 
be used for information processing. The downside of it is 
that it is in principle impossible to have maximally en- 
tangled states, and consequently there is an additional 
source of imperfections. Furthermore, finding concentra- 
tion schemes in the continuous variable regime is made 
particularly difficult by a no-go theorem that states that 
distillation of Gaussian states cannot be carried out us- 
ing Gaussian operations alone [4-6 . One has to resort 
to challenging non-Gaussian operations, like strong high- 
order nonlinear interactions or photon counting. 

It is worth noting that in the continuous variable 
regime, more or less all of the proposed schemes for en- 
tanglement distillation contain variations of a common 
theme: beamsplitter interactions. On the theory side, 
Opartny et al. ^ showed how beamsplitters and pho- 
ton subtraction can be used to increase the entangle- 
ment in a two mode squeezed vacuum (TMSV) as part 
of a teleportation scheme. The distillation aspect of this 



was seized upon and improved [SHTTj and a rigorous the- 
oretical description was given by Kitagawa et al. |12j . 
Experimentally, there have been some notable successes 
[131 [13], including with non-Gaussian noise [T^. In a 
fairly recent proposal [16 a complete GV entanglement 
"distillery" has been proposed, harnessing limited phys- 
ical space for storing quantum states and distilling en- 
tanglement. There, clever manipulation of an imperfect 
quantum memory complements the beamsplitter based 
entanglement concentration scheme. 

However, the net could be cast wider. Theoretical 
work has been carried out utilizing other interactions be- 
tween an entangled state and ancillary modes, notably 
the cross-Kerr effect [T71 [TH] , in which the nonlinear in- 
teraction plays the role of the necessary non-Gaussian 
element, and weak values [19]. Still, very little is known 
about entanglement distillation protocols using so-called 
Quantum Non-Demolition (QND) interaction Hamilto- 
nians, although they lend themselves readily to many 
experimental systems where distillation could be useful. 
For example the off-resonant dipole interaction between 
heavily polarized light and the macroscopic spin proper- 
ties of an atomic ensemble can be modeled as QND and 
two atomic ensembles were successfully entangled, albeit 
weakly, in this way [SHj . Furthermore, for light it was also 
shown that it is possible to create a QND Hamiltonian 
using optical instruments such as biased beamsplitters 
and squeezers in combination [21j . 

The Quantum Non-Demolition interactions |22', '23' are 
well known in quantum optics for quantum metrology, 
entanglement generation and other purposes. They are 
described by Hamiltonians of the form e.g. 



kPsPa 



(1) 
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where k is the interaction strength and Ps/a a-re the mo- 
mentum quadratures of the system and ancillary mode. 
The system and ancillary mode interact in such a way 
as to leave one quadrature component of each subsystem 
intact whilst phase shifting the conjugate components. 

In this paper, we suggest and explore two procrustean 
entanglement distillation schemes that use QND interac- 
tions, and quantify their performance showing the (prob- 



2 



abilistic) increase of the the logarithmic negativity of a 
two mode squeezed vacuum, defined in the Fock basis as 



(2) 



where |n) is the photon number and A = tanhr is the 
squeezing. For the quantitative analysis of these schemes, 
we amalgamate ideas on using the well-known formal- 
ism of Gaussian states to represent non-Gaussian states 
[Mj (used more recently in [5S]) and calculate density 
matrix elements using multivariate Hermite polynomials 
|26H28j ■ Our approach is incredibly versatile, making ad- 
ditional interactions and double passes of the interacting 
modes as simple as multiplying 8x8 matrices and al- 
lowing for detector inefficiencies to be considered simply. 
Hence numerical calculations are made fast and accurate. 
Thus the novelty of the proposed approach to distillation 
lies both in the scheme design (QND-based) and in the 
methodology of the calculations. Notably, the suggested 
scheme is not bound to a certain physical system and 
can be applied to the distillation of entanglement in light 
modes, atomic ensembles or in hybrid systems, the essen- 
tial requirement being solely the form of the Hamiltonian 
underlying the interaction. 

The paper is organized as follows. In section [Tl] we 
remind the reader of the important properties of Gaus- 
sian states required. In section III we study the first 



scheme, a procrustean entanglement concentration pro- 
tocol replacing the conventional beamsplitter interaction 
with the QND interaction between a Gaussian ancilla 
and a Gaussian entangled state. More precisely, we cal- 
culate the logarithmic negativity of the non-Gaussian 
mixed state resulting from the QND interaction between 
the two mode squeezed vacuum and a Gaussian ancil- 
lary mode subsequently subjected to a non-Gaussian 
measurement. In section |IV[ we explore an alternative 
scheme where the ancillary mode is non-Gaussian and the 
post-measurement is Gaussian. In both cases the entan- 
glement of the system is shown to increase demonstrating 
a successful continuous variable entanglement distillation 
based on the QND Hamiltonian. 



II. PRELIMINARIES 

An N-mode Gaussian state can be completely defined 
by the first and second order moments of the quadra- 
ture operators Rj with R = . . . ,xn,Pn)'^ ■ The 
Wigner function for such a state is given by 



W(a:i,Pi, . . .,Xm,Pn) = 



exp 



R^_d^j 7-1 (R-d) 
(3) 

where d is the displacement vector d — 
{{xi) , {pi) , . . . , {xjy) , {pn))'^ and 7 is the covariance 



matrix of the state in question, defined as 

llm — ( RlRm + R-mRl ) ~ '^didm- 



(4) 



The two mode squeezed state that we wish to distil has 
covariance matrix 



TMSS 

Jab 



I cosh(2r) sinh(2r) 

cosh(2r) - sinh(2r) 

sinh(2r) cosh(2r) 

V -sinh(2r) cosh(2r) 



. . (5) 

This Gaussian distribution is centred at the origin in 
phase space and so d of equation ([s]) is the zero vector. 
The squeezing parameter is denoted r. 

As a well-known consequence of the Stone- Von- 
Neumann theorem, any quadratic unitary operation at 
the Hilbert space level corresponds to a symplectic op- 
eration at the phase space level, that preserves the sym- 
plectic form — (SjLiiay. In what follows, the following 
symplectic operations are relevant. Firstly, the QND in- 
teraction between two modes (system and ancilla) can be 
described by its action on phase space, e.g. 



S, 



(nPP) 
QND 
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(6) 



corresponds to the Hamiltonian of the form given in 



Eq. (1). Similarly, S, 



a, 



QND 



QND 



and Sgj^f^^ can 



also be defined where X is the position quadrature. The 
phase shift of one mode is defined as 



SPH (0) = 



cos 6 sin 9 
— sin 9 cos 9 



(7) 



where 9 is the angle through which the Gaussian mode 
is rotated in phase space. We also define the single mode 
squeezer 



SsQ I 







(8) 



which yields jsQ = diag (e^*,e~^*) when applied to the 
vacuum covariance matrix ^yac = 1 . If s > then the 
mode is squeezed in p and the mode is squeezed in x if 
s < 0. 



Noisy quantum channels and phase insensitive ampli- 
fiers can not be modeled by Gaussian unitary transforma- 
tions. One of the biggest losses by detectors is when they 
erroneously fail to register a detection. For this reason we 
model inefficient detectors as being perfect detectors be- 
hind a beamsplitter of transmittance t] (see e.g. [9|, I29j). 
This can be seen as combining the mode with the vacuum 
on a beamsplitter and tracing out the vacuum. Such an 
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Figure 1: The two modes of a two mode squeezed vacuum in- 
teract with two vacuum modes which are subsequently mea- 
sured for the presence or absence of photons. If successful, 
the entanglement of the two mode squeezed state increases. 

operation transforms the covariance matrix as 

7^5io3.(»7)7Cs(^) + (l-'?)l (9) 
where ^loss iv) = V^'^- 



That is, the total Wigner function Wab {xa,Pa,xb,Pb) 
is given by a superposition of normalized Gaussian 

Wigner functions W^';^'' = V^ab (^AiPa,xb,Pb) with 
i = meaning that the mode a has been projected onto 
|0) (0| and i — 1 imphes that a has been projected onto 1 
(i.e. traced out). The index j tells what has happened to 
mode b. The constants p(*J) are the probabilities associ- 
ated with W^'^ and M is a global normalisation constant 

given as M = {J2^,J (-1)'^' P^"'^) 

The 4x4 covariance matrices j\g that make up the 
Gaussian Wigner functions yVAB ^^"^ their normalisa- 
tion constants are calculated from 7 by consider- 
ing the Wigner Overlap formula. Tracing out a mode 
corresponds to an overlap with Wi = (l/27r) and a pro- 
jection onto the vacuum corresponds to an overlap with 
Wo = TT^^ exp[— a;^ — P^]- For further details see e.g. [24] . 



We define 



III. PROTOCOL I: GAUSSIAN ANCILLARY 
MODE AND NON-GAUSSIAN MEASUREMENT 



The first protocol, depicted in Figure [T] facilitates the 
concentration of the entanglement of some continuous de- 
gree of freedom coupled to ancillary modes with QND in- 
teractions. The entangled state is in the form ([2|, where 
the exact meaning of the basis states \nn) is determined 
by the used system. The ancillary modes are in a vacuum 
state and so are described in phase space by a Gaussian 
Wigner function. The required non-Gaussian element |4]- 
|6] is in this case fulfilled by a non-Gaussian measure- 
ment. Thus we begin with a two mode squeezed vacuum 
(modes A and B) and two vacuum ancillary modes (a 
and b) (Figure Jl|. The initial covariance matrix is given 



by 7init = 1a^^ ® lvac,a © Tvac,b- After the desired 
QND interaction is performed between A and a, and B 
and b, the two ancillary modes are detected for the pres- 
ence or absence of photons on a detector of efficiency 77. 
That is, the non-Gaussian measurement 11 = 1 — |0)(0| is 
performed on modes a and b. Immediately prior to the 
measurement, the covariance matrix of all four modes 
can be described by 

7 = 

{IaB © Sloss,ab (j?)) S'qND 7init -Sq^d i'^AB © S'loss.ab {v)) 

+ (C ® (1 - 77) lafc) . (10) 

Conceivably, other interactions could be considered by 



modifying (10), such as double-pass schemes (see [TTl 
l30]). After the measurement, the total Wigner function 
of modes A and B is given by 



Wai 



M 



1 

i,j=0 



(11) 



7 



Tab o- 



(12) 



where Tab describes the entangled modes A and B and 
Tab describes the ancilla modes a and b with a captur- 
ing the cross correlations between entangled and ancilla 



modes. Then j^^b 



^ AB) 



1 AB - 1 



where 



ab 



<i3) 
ab 



Tab + il-i)l®{l~j)l. 



(13) 



The normalisation constants P^'J) can be calculated by 

(ii) 

the rather straightforward integration of W\g which can 
simply be thought of as the overlap of modes A and B 
with Wi. Consequently, 



22-*-Vdetr 



detr^i]Jdetri-) 



(14) 



At this point, we have specified all the quantities in 
the expression (11) for the Wigner function of the two 



mode squeezed state after QND interactions with vac- 
uum modes that are later detected for the presence or 
absence of photons on inefficient detectors. Hence we are 
able to describe the entangled state modified by the QND 
interaction by use of the Wigner function ( 11 ). However, 



to show that the entanglement in the two modes is in- 
creased, we require the density matrix elements of the 
state. These can be calculated with the help of the Q 
function and multidimensional Hermite polynomials |26l - 



The Q function is defined as the convolution of the 
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Wigner function with the vacuum [3T] and so 



q(u) 



exp 



7e 



TT^ydct (7^) + 1) 



(15) 



where 7?. = (a, a^, /3, and U is the unitary transfor- 
mation between {xa,Pa,xb,Pb)^ (^a,a\b,b^^ : 



U = 




(16) 



The Q function then acts as a generating function for the 



density matrix elements, pl 



(ii) 



(2 



(;| (m| p^») \n) \p) 
Q(j^i(a,/3)eH'+l^l " 



la=/3=0 



(17) 



The matrix element is most easily expressed in terms 
of the four dimensional Hermite polynomials Himnp (see 
appendix []) . After successful detection the total density 
matrix is given by p with matrix elements 



4M 



Plmnp 



1 



Ininp 



^l\m\n\p\ f-^^ ^det {^'^^1) + l) 



(18) 



where for brevity we have defined ii\mnp '"^^ = 

^/mrir °^ (0, 0, 0, 0). The matrices C^*-'^ are derived from 
the matrix equation 



C^^^) = B 



1 



D 



(19) 



with 
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1 





/ 



(20) 



The purpose of B and D is simply to rearrange the ele- 
ments of the matrix in such a way as to make the matrix 
compatible with the form given in (A.l) and to make 

Plmnp proportional to i?;„„p . For single QND passes 
the matrix C'-*-''' is always real. 



With the density matrix elements of modes A and B 
we can now show that the entanglement of the system 
has increased. To quantify this we use the additive Log- 



Figure 2: Logarithmic Negativity of TMSV (solid line), and 
after a QND interaction of the form kXP with strength 0.1 
(dotted), 0.5 (dashed), and 1 (thick dashed). The scheme used 
is that depicted in Figure [T] The entanglement increase is 
best when weak interactions are used. Strong interactions can 
decrease the entanglement in the initial TMSV. The vertical 
lines from left to right indicate where the initial squeezing r 
is IdB, 3dB, 5dB, and 7dB. 



arithmic Negativity [32] defined as 

£{pab)^\^\p'ab \ -ln|2AA(p^B) + l| 



(21) 



where ||/5}i^|| trace- norm of the partial transpose 

of pAB- The Negativity, M is defined as the sum of the 
negative eigenvalues of 



\Pab\ 



The logarithmic negativity of the mixed state result- 
ing from this protocol is shown in Figure [2l The graph 
demonstrates the performance of the protocol in different 
regimes: depending on the initial entanglement content 
in the two-mode squeezed vacuum and the QND inter- 
action strength (see caption). Successful entanglement 
concentration is achieved for weak to medium values of 
K, i.e. for a moderate QND coupling between the TMSV 
and ancilla modes, for a wide range of the initial squeez- 
ing in the TMSV. If the interaction is too strong (e.g. 
K = 1) then the effect of the QND interaction on the 
entanglement content is negative, but as the interaction 
strength decreases, the logarithmic negativity of the re- 
sultant state increases and the performed QND operation 
is beneficial at least when the initial squeezing of the 
TMSV is weak. That is, as the interaction strength de- 
creases, it becomes more likely that just a single photon 
is subtracted from each mode of the two mode squeezed 
vacuum and these photons are then detected. As ex- 
pected, this compares well to the case where beamsplit- 
ters are used in place of QND interactions, investigated 
numerically in [II] for light TMSV and in [11^ for TMSV- 
entangled atomic ensembles. After all, if the pure state 
resulting from a single photon subtraction in both modes 
of the entangled state is more entangled than the initial 
state, then it should not matter which interaction exactly 
is used to subtract the photons. 
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IV. PROTOCOL II: NON-GAUSSIAN 
ANCILLARY MODE AND GAUSSIAN 
MEASUREMENT 

Whereas Protocol I relies on a Gaussian ancillary mode 
interacting with the Gaussian TMSV via QND interac- 
tion, and a subsequent non-Gaussian measurement, we 
here consider a protocol which uses a non-Gaussian an- 
cilla and a Gaussian measurement after the QND inter- 
action (Figure [3| . The advantage of this new protocol is 
that probabilistic non-Gaussian measurements are per- 
formed in the preparatory stage and, if successful, the 
QND interaction can be carried out in the next step. 
Also, the detection part of the scheme outperforms the 
first protocol as the Gaussian measurement, homodyne 
detection, can be highly efficient. The same methods as 
above can be used to calculate the logarithmic negativity 
of this scheme. 



A. Preparatory step 

In the protocol of Figure [3] the ancillary modes a and 
h are in a quantum state with a non-Gaussian Wigner 
function. Specifically, they are in a photon-subtracted 
squeezed vacuum state. We consider first the prepara- 
tion of the ancillary modes. In contrast to the previous 
scheme the modes a and h are initially in a squeezed 
vacuum state, ^ab = lSQ,a (s) ® 7SQ,b (s)- The removal 
of photons from the mode is modeled in the following 
way. The modes a and h pass through a beamsplitter of 
transmittance and any subtracted photons are subse- 
quently detected with efficiency rj. As the transmittance 
increases, the model approaches the limit of just a 
single photon being subtracted from the mode. In what 
follows, in all numerical calculations we take T — 0.95. 

If photon subtraction is successful, the resultant 
non-Gaussian Wigner function becomes a superposi- 
tion of Gaussian Wigner functions in the form Wat = 
MabY.^,J{-^Y^' P^^'^y^ai^ with probability Mab = 
(— 1)*^"' . The corresponding covariance ma- 
trices that make up the Wigner functions W^f ■* are given 

by lab = 7prep ® 7prep whcrC 



/ prep 







(22) 



and 



(2 - k) (e 



±2s 



1 



1 + (1 - fc) (1 + (e±2« - 1) (1 - ^2)) • 



The coefficients are calculated as 

22— J- V^+(*)^_(*)^+(j)^_(j) 



^/T+{i)T-{i)T+{j)T^{j) 



(23) 



(24) 



with 



T±{k) = \-k- 



1 + (e±2« - 1) t2 



(1 - ry) (e 



±2s 



1) (1 -T2) 



,±2s 



(25) 



At this point we have a non-Gaussian state Wab, consist- 
ing of a superposition of Gaussian states. If the prepara- 
tory step could be performed independently and only 
passed to the interaction and homodyne measurement 
if successful then we could consider -q = 1. 



B. Interaction and Detection 

A non-Gaussian state produced in the preparatory step 
further interacts with the TMSV via QND coupling and 
the ancilla is subsequently passed on for homodyning on a 
detector of efficiency ^. An angle 6 describes the phase of 
the homodyne measurement, that is the angle at which 
homodyning is performed in phase space. An angle of 
= corresponds to a measurement of position quadra- 
ture xg^Q = X (the X marginal distribution) whereas an 
angle of 7r/2 describes the momentum quadrature mea- 
surement Xg^^f2 = P (the p-distribution). Thus the ho- 
modyne measurement is characterised by the generalized 
quadratures xe^a,xe,b- Directly before measurement, we 
can calculate the covariance matrices to be 



lABab 



with 



{Iab © SpH {9)) (1 © a) (Iab ® {9)) 



(26) 



Tint - JQND 



TMSS ^ Jij) 

Iab © lab 



"^QND- 



(27) 



To see what happens after the homodyne detection we 
project the quadratures xe^a and xe^b onto an outcome z. 
We also trace out the conjugate quadratures {xg^^/2.a 
and Xg+^/2,b)- 

At this point the remaining covariance matrix for A 
and B and the correlations due to the z measurements 
can be described by a 6 x 6 matrix such that 



(28) 



where .4^'-'^ is a 4 x 4 matrix (modes A and B), S^'-J) is a 
2x2 matrix descibing the z correlations, and C^*^-* con- 
tains the cross correlations. The probability of projection 
of modes a and b is given by 



X exp 



(z,z) 



(z,z) 



(29) 
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Figure 3: A non-Gaussian ancillary state is first created in a preparatory step by subtracting photons from a squeezed vacuum. 
The ancilla then interacts with the TMSV via QND interactions and is measured by homodyne detection. 
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Figure 4: Performance of protocol II (Fig. |3|: Logarithmic 
Negativity of TMSV before (solid line), and after a QND in- 
teraction for different initial squeezing of the ancilla modes. 
The QND interaction is of the form kXP and has strength 
K = 0.5; homodyne outcome z « 0. The ancillary modes are 
squeezed by s = 0.1(0.87dB) (black dotted), s = 0.5(4.34dB) 
(black dashed) and s = l(8.69d_B) (thick black dashed). The 
improvements of entanglement in the TMSV gained by in- 
creasing the squeezing of the ancillary modes is small when 
the initial squeezing of the TMSV r is small. The vertical 
lines from left to right indicate where the initial squeezing of 
the TMSV r is IdB, 3dB, 5dB, and 7dB. 77 = ^ = 1. 



Figure 5: Performance of protocol II (Fig. |3|: Logarithmic 
Negativity of TMSV before (sohd line), and after a QND in- 
teraction for different QND interaction strengths. The QND 
interaction is of the form kXP. The non-Gaussian ancillary 
modes are initially squeezed by s = 0.2(1.74ci_B) and it is as- 
sumed that the homodyne measurement outcome is 2: ~ 0. 
The interaction strengths are k = 0.1 (black dotted), k = 0.5 
(black dashed) and k = 1 (thick black dashed). An intermedi- 
ate interaction strength is preferable. The vertical lines from 
left to right indicate where the initial squeezing of the TMSV 
r is IdB, 3dB, 5dB, and 7dB. ri = ^ ^ I. 



and the Q function describing modes A and B reads 



Vdet -A('j),t,('J)-iA('^)t 



X exp 



where 



(30) 



(31) 



A('J') = 2 (z, z) Cfe)'^yt(*^')-i (yi('^')-i +1) ^ (32) 

By defining = B - (1/2) l] D and A' = AD 

we can write the density matrix elements as 



a// m T7, n — ' 



Vde^e ^ i?/*„p -'^ > 

(33) 
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with 



DISCUSSION AND CONCLUSIONS 



t 



(34) 



From this, the logarithmic negativity of modes A and B 
can be calculated as before. 

In contrast to the first protocol, the second protocol is 
sensitive to which interactions are used. The position and 
momentum correlations in the TMSV are not mixed by 
the QND interactions and so the effects induced by inter- 
action Hamiltonians e.g. h[^^^^ and h[^^^^ are equiva- 
lent. Correspondingly, the choice of measurement on the 
ancillary modes is important. If, for example, the inter- 
action H^^^^^ is used, then the momentum quadratures 
of the ancillary modes are, by definition, unaffected but 
information about the TMSV is imprinted on the posi- 
tion quadrature distribution. If homodyne measurements 
on the ancillas are performed owp {6 = 7r/2) then the re- 
sulting outcome can reveal nothing about the state of the 



TMSV and the probabilities (29) are independent 
of K and r. As nothing can be learned probabilistically 
about modes A and _B, procrustean entanglement con- 
centration cannot occur. The effects on the TMSV can 
be transformed away by local Gaussian operations and 
ancillary modes and so the entanglement of the TMSV 
is unchanged. The above is true irrespective of whether 
the ancillary modes are initially squeezed in position or 
momentum. 

If homodyne measurements were taken of the position 
quadratures {9 — 0) then information about the TMSV 
will have been probabilistically imprinted on this mea- 
sured distribution. If h[^^^^ is used then the momen- 
tum quadratures of the TMSV contain information about 
the momentum quadratures of the ancillary modes. That 
is, some noise has been added to modes A and B which 
can assist or disrupt entanglement concentration. If the 
ancillary modes are squeezed in momentum, then only a 
little noise is added to the p quadratures of the TMSV. 
The result is an increase in entanglement dependent on 
s and K (Figure |4| . If the ancillary modes are squeezed 
in position then the momentum quadratures are anti- 
squeezed and so a lot of noise is added to the TMSV, 
having a detrimental effect on the entanglement. 

As can be seen in Figures |4] and [5] successful entangle- 
ment concentration is achieved if the ancillary modes are 
squeezed in p-quadrature, the QND interaction is of the 
form H^^^'^^ , and homodyne measurement is performed 
in x-quadratures. The protocol is noticeably insensitive 
to ancillary mode squeezing - the logarithmic negativity 
of the resulting state is largely unaffected for low levels 
of TMSV initial squeezing. 

Unlike Protocol I a weak interaction strength ruins the 
entanglement in the system. Similarly, if the strength is 
too strong then the entanglement decreases. If k w 0.5 
then concentration successfully occurs. 



We have presented here two procrustean entangle- 
ment concentrations schemes utilizing Quantum non- 
Demolition (QND) interactions and photon detectors. 
The first scheme relied upon QND interactions between 
Gaussian ancillary modes and a TMSV to successfully 
subtract photons from the TMSV as heralded by on/off 
detectors. We have shown how to efficiently calculate 
the density matrix elements of the resulting quantum 
state, which can then be used to calculate the logarith- 
mic negativity of the state. This is a non-trivial task. 
In the asymptotic limit of the QND interaction strength 
K — 0, one would find that a single photon is subtracted 
from each arm of the TMSV, although the probability of 
heralding entanglement concentration approaches zero. 
This is intuitively correct as the behavior mimics the 
protocol of [H] in which highly transmissive beamsplit- 
ters are used in place of weak QND interactions. There, 
as the transmittance approaches 1, one finds that a single 
photon is subtracted from each arm, which is the optimal 
outcome for entanglement concentration. The results of 
Browne et al. can be revisited by replacing Sqnd with 
•^BS (T) in equation ( 10 ) where Sbs (T) is the symplectic 
operation corresponding to a beamsplitter transforma- 
tion. For two modes this is 



Sbs (T) 



It Vi -^2 \ 

T V1-T2 

^/\^^ -T 

V Vi-r2 -T ) 

(35) 

We find that the method we have used for calculating 
density matrix elements for this protocol offers an im- 
provement in numerical speed and efficiency over calcu- 
lations in the Fock basis directly. 

Double pass schemes, i.e., letting the ancillae interact 
with the TMSV twice, do not give any advantage in pro- 



tocol I. For example s'^^^^^ s'^^^^ 
the logarithmic negativity is highest for weak interaction 
strength. By altering the interaction between passes e.g. 

'^qnd'^^'^qnd^^ there is some advantage over a double 
pass scheme using the same interaction twice but this is 
still eclipsed by the single pass schemes. 

As to detector efficiencies, at low levels of initial TMSV 
squeezing r, there can still be an increase in entanglement 
for weak interactions such as k = 0.1 when the detector 
efficiency exceeds approximately 50%, r/ > 0.5. 

The second protocol relied on QND interactions to 
mix a TMSV with photon subtracted squeezed vacuum 
modes, that is with a non-Gaussian ancillae. The non- 
Gaussian ancillary modes are then detected after the in- 
teraction using homodyne detectors. As was stated in 
the previous section, the success of this scheme is de- 
pendent on the initial squeezing of the ancillary modes, 
the choice of interaction, and the angle Q of homodyning. 
It is most likely that a homodyne measurement yields a 



_ c.((ki-I-K2)XP) , 

— '-'qnd 
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result z = 0, and the logarithmic negativity shows an 
improvement so long as the measurement outcome does 
not stray too far from this. 

The benefit of the second protocol is that the prob- 
abilistic non-Gaussian state preparation could be per- 
formed off-line and therefore the efficiency 77 of the on/off 
detectors could be assumed to be good. Then only the 
homodyne detector efficiency needs to be taken into ac- 
count. 

Let us now consider the feasibility of the implemen- 
tation of the protocols presented here in a laboratory 
setting. In a purely optical setup in which the TMSV is 
created by parametric down conversion and the QND op- 
erations are performed with beamsplitters and squeezers, 
protocol I is probably the easiest to implement, although 
other effects such as a nonzero dark count rate of the 
detectors would need consideration. 

However, in the atomic systems of [20] where the en- 
tangled macroscopic spin states of two cesium gas sam- 
ples represent the TMSV and off-resonant dipole interac- 
tions with strongly polarized light form the QND inter- 
actions, the vacuum ancillary modes are not true vacuum 
modes. They are instead modes of heavily polarized light. 
The heralding of the entanglement concentration comes 
from detecting photons for which the interaction has al- 
tered the polarization, but this requires heavy filtering, a 
problem that is most likely insurmountable with current 



technology. Protocol II does not suffer this problem and 
so is the best choice in the atomic case. The probabili- 
ties of success in both cases are comparable with those 
for the beamsplitter-based light schemes, already demon- 
strated [13J in the laboratory. The possibility of distilling 
entanglement in atomic ensembles represents the main 
motivation for using a general QND Hamiltonian for en- 
tanglement concentration. 
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Appendix: Multivariate Hermite Polynomials 

In order to calculate the density matrix elements of 
a Gaussian state, the multivariable Hermite polynomials 
[33| are an invaluable tool. The four dimensional Hermite 
polynomials are defined by 



H. 



{e,A} 



(yi,2/2,y3,y4) 



\r~\-s~\-t-\-v 



exp [y'^ey + Ay] 



gr QS Qt QV 



where y = (2/1,2/2,2/3,2/4) , O is a 4 x 4 matrix and 
A = (Ai, A2, A3, A4). The matrix 9 is related to the 
inverse of the covariance matrix of the Gaussian state 
and the matrix A is related to any displacements. The 
density matrix elements are found by evaluating the Her- 
mite polynomials at yj = 0. A recursion formula can 
be derived to help speed up the numerical calculation of 
the multivariable Hermite polynomials by directly using 
( A.l) and substituting there e.g. r + 1. This gives 



{0,A} 



r-\-l,s,t,v 



{e,A} 

r,s,t,v 



t (613 + 631) ^f;,'^^,^ 



(A.2) 



{e,A} _ „{e,A} 



TJ; , (0, 0, 0, 0). Similar for- 



In the above -ff, 

mulae can be derived for Hj:^^\ „ etc. by replacing the 



coefficients in (A.2 1. There are some symmetries which 



can be exploited. As the resulting density matrix is nec- 
essarily Hermitian, Hj:^'^^ = Hf\fl'^\ Furthermore, if 
A = as in protocol I or in protocol II if the homodyne 
projection yields z — 0, then only Hermite polynomials 
for which the sum of the indices is an even number need 
to be calculated as whenever the sum is odd, the corre- 
sponding polynomial is zero. 
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